Figure1_S1-S6

Tom LaSalle

This document contains all the code necessary to generate the plots for Figure 1 and related supplementary figures (S1-S6). Plots are subsequently edited in Adobe Illustrator to produce the final figures.

Load the necessary libraries:

library(knitr)
library(ggplot2)
library(ggrepel)
library(RColorBrewer)
library(plyr)
library(dplyr)
library(DESeq2)
library(openxlsx)
library(corrplot)
library(ggmosaic)
library(cowplot)
library(Seurat)
library(heatmap3)
library(ggpubr)
library(umap)
library(fgsea)

Import the metadata and keep only data for which neutrophil enrichment was performed:

prefix <- "~/Downloads/Github/"
metadata_bypatient <- read.xlsx(paste0(prefix,"Tables/TableS1.xlsx"), sheet = 2)
qc_data <- read.xlsx(paste0(prefix,"Tables/TableS1.xlsx"), sheet = 9)
metadata_bypatient <- metadata_bypatient[which(metadata_bypatient$Public.ID %in% qc_data$Public.ID),]
genomic_signatures <- read.xlsx(paste0(prefix,"Tables/TableS1.xlsx"), sheet = 10)
metadata_long <- read.xlsx(paste0(prefix,"Tables/TableS1.xlsx"), sheet = 4)

Due to an institutional IRB-approved waiver of informed consent, all clinical data are reported in quintiles, and not all clinical characteristics are publicly available. Correlations were calculated between the following clinical variables: Age, Sex, BMI, Heart condition, Lung condition, Kidney condition, Diabetes, Hypertension, Immunocompromised, Symptom duration, Respiratory symptoms, Fever symptoms, GI Symptoms, Acuity (D0, D3, D7, Max within D7-D28, D28, Max Overall), ANC (D0, D3, D7), ALC (D0, D3, D7), AMC (D0, D3, D7), Creatinine (D0, D3, D7), CRP (D0, D3, D7), D-dimer (D0, D3, D7), LDH (D0, D3, D7), Troponin detected within 72 hours, Intubation status, Race (White, Black, Asian, Other), Hispanic ethnicity, CXR infiltrates, and Smoking history. Multiple hypothesis correction was performed to determine significant correlations, and then those parameters with significant correlations were plotted in a correlation heatmap. Though it is not possible to show all of the data, the code to generate the figure is below and is written as if the additional columns were present in the metadata file:

# # Reverse factor levels of Acuity such that positive correlations with Acuity indicate correlations with worse disease severity
# metadata_bypatient$Acuity.0 <- mapvalues(metadata_bypatient$Acuity.0, from = c(1,2,3,4,5), to = c(5,4,3,2,1))
# metadata_bypatient$Acuity.3 <- mapvalues(metadata_bypatient$Acuity.3, from = c(1,2,3,4,5), to = c(5,4,3,2,1))
# metadata_bypatient$Acuity.7 <- mapvalues(metadata_bypatient$Acuity.7, from = c(1,2,3,4,5), to = c(5,4,3,2,1))
# metadata_bypatient$Acuity.28 <- mapvalues(metadata_bypatient$Acuity.28, from = c(1,2,3,4,5), to = c(5,4,3,2,1))
# metadata_bypatient$Acuity.7.to.28 <- mapvalues(metadata_bypatient$Acuity.7.to.28, from = c(1,2,3,4,5), to = c(5,4,3,2,1))
# metadata_bypatient$Acuity.max <- mapvalues(metadata_bypatient$Acuity.max, from = c(1,2,3,4,5), to = c(5,4,3,2,1))
# 
# rownames(metadata_bypatient) <- metadata_bypatient$Public.ID
# 
# metadata_bypatient <- metadata_bypatient[,-which(colnames(metadata_bypatient) %in% c("Study_ID","Public.ID","D0_draw","D3_draw","D7_draw","DE_draw","D0_Neu_Sample","D3_Neu_Sample","D7_Neu_Sample","DE_Neu_Sample"))]
# metadata_bypatient <- metadata_bypatient[metadata_bypatient$COVID == "1",]
# metadata_bypatient <- metadata_bypatient[,-which(colnames(metadata_bypatient) %in% c("COVID"))]
# 
# cormatrix <- matrix(0L, nrow = ncol(metadata_bypatient), ncol = ncol(metadata_bypatient))
# pmatrix <- matrix(0L, nrow = ncol(metadata_bypatient), ncol = ncol(metadata_bypatient))
# rownames(cormatrix) <- rownames(pmatrix) <- colnames(cormatrix) <- colnames(pmatrix) <- colnames(metadata_bypatient)
# 
# for (i in 1:nrow(cormatrix)){
#   for (j in 1:nrow(cormatrix)){
#     stats <- cor.test(x = as.numeric(metadata_bypatient[,i]), y = as.numeric(metadata_bypatient[,j]), use = "pairwise.complete.obs", method = "kendall")
#     cormatrix[i,j] <- stats$estimate
#     pmatrix[i,j] <- stats$p.value
#   }
# }

cormatrix <- as.matrix(read.xlsx(paste0(prefix,"Tables/TableS1.xlsx"), sheet = 5, rowNames = TRUE))
pmatrix <- as.matrix(read.xlsx(paste0(prefix,"Tables/TableS1.xlsx"), sheet = 6, rowNames = TRUE))

ptable <- matrix(0L, nrow = nrow(pmatrix)^2, ncol = 3)
colnames(ptable) <- c("Row","Column","pvalue")
ptable[,1] <- rep(colnames(pmatrix), each = nrow(pmatrix))
ptable[,2] <- rep(colnames(pmatrix), nrow(pmatrix))
for (i in 1:nrow(ptable)){
  ptable[i,3] <- pmatrix[rownames(pmatrix) %in% ptable[i,1],colnames(pmatrix) %in% ptable[i,2]]
}
ptable <- as.data.frame(ptable[complete.cases(ptable[,3]),])
ptable$fdr <- p.adjust(ptable[,3], method = "fdr")
ptable <- ptable[ptable[,4] < 0.05,]
ptable <- ptable[ptable[,1] %in% c("abs_neut_0_cat","abs_neut_3_cat","abs_neut_7_cat") | ptable[,2] %in% c("abs_neut_0_cat","abs_neut_3_cat","abs_neut_7_cat"),]

cormatrix <- cormatrix[rownames(cormatrix) %in% unique(c(ptable[,1],ptable[,2])),colnames(cormatrix) %in% unique(c(ptable[,1],ptable[,2]))]
pmatrix <- pmatrix[rownames(pmatrix) %in% unique(c(ptable[,1],ptable[,2])),colnames(pmatrix) %in% unique(c(ptable[,1],ptable[,2]))]

rownames(cormatrix) <- colnames(cormatrix) <- rownames(pmatrix) <- colnames(pmatrix) <- mapvalues(rownames(cormatrix), from = c("BMI.cat","abs_mono_0_cat","abs_mono_3_cat","abs_mono_7_cat","Age.cat","creat_0_cat","creat_3_cat","creat_7_cat","ldh_0_cat","ldh_3_cat","ldh_7_cat","crp_0_cat","crp_3_cat","crp_7_cat","Acuity.28","Acuity.7","Acuity.3","Acuity.0","Acuity.max","abs_neut_0_cat","abs_neut_3_cat","abs_neut_7_cat","ddimer_0_cat","ddimer_3_cat","ddimer_7_cat"), to = c("BMI","AMC, D0", "AMC, D3","AMC, D7","Age","Creatinine, D0","Creatinine, D3","Creatinine, D7","LDH, D0","LDH, D3","LDH, D7","CRP, D0","CRP, D3","CRP, D7","Acuity, D28","Acuity, D7","Acuity, D3","Acuity, D0","Acuity, Max","ANC, D0","ANC, D3","ANC, D7","D-dimer, D0","D-dimer, D3","D-dimer, D7"))

Figure 1B:

corrplot(cormatrix, method = "square", type = "lower", order = "hclust", hclust.method = "ward.D", tl.col="black", col=colorRampPalette(c("blue","white","red"))(200), p.mat = pmatrix, sig.level = 0.05, insig = "blank")

Next we highlight the correlations between ANC on each individual day and Acuity Max, using Kendall’s tau correlation:

metadata_long <- metadata_long[which(metadata_long$Public.ID %in% qc_data$Public.ID),]
metadata_temp <- metadata_long[complete.cases(metadata_long$ANC.matched) & metadata_long$COVID == 1,]

metadata_temp$Acuity.max <- factor(metadata_temp$Acuity.max, levels = c("5","4","3","2","1"))
corvals0 <- cor.test(x = as.numeric(metadata_temp[metadata_temp$Day == "D0",]$Acuity.max), y = as.numeric(metadata_temp[metadata_temp$Day == "D0",]$ANC.matched), use = "pairwise.complete.obs", method = "kendall")
corvals3 <- cor.test(x = as.numeric(metadata_temp[metadata_temp$Day == "D3",]$Acuity.max), y = as.numeric(metadata_temp[metadata_temp$Day == "D3",]$ANC.matched), use = "pairwise.complete.obs", method = "kendall")
corvals7 <- cor.test(x = as.numeric(metadata_temp[metadata_temp$Day == "D7",]$Acuity.max), y = as.numeric(metadata_temp[metadata_temp$Day == "D7",]$ANC.matched), use = "pairwise.complete.obs", method = "kendall")

metadata_temp$ANC.matched <- factor(metadata_temp$ANC.matched, levels = c("5","4","3","2","1"))

my.cols <- brewer.pal(9, "RdBu")
p1 <- ggplot(data = metadata_temp[metadata_temp$Day == "D0",]) + geom_mosaic(aes(x = product(ANC.matched, Acuity.max), fill = ANC.matched), offset = 0.005) + coord_fixed(ratio = 1) + theme_bw() + scale_fill_manual(values = (my.cols[(c(1,2,3,7,9))])) + ggtitle("Day 0") + annotate(geom="text", x=0.5, y=.8, label=paste0("τ =",round(cormatrix[rownames(cormatrix) == "ANC, D0",colnames(cormatrix) == "Acuity, Max"],3),"\npadj =",round(pmatrix[rownames(pmatrix) == "ANC, D0",colnames(pmatrix) == "Acuity, Max"],12)), color="black") + theme(legend.key.size=unit(4,'mm'), legend.text=element_text(size=5), legend.title=element_text(size=5), axis.title = element_text(size = 7))
p2 <- ggplot(data = metadata_temp[metadata_temp$Day == "D3",]) + geom_mosaic(aes(x = product(ANC.matched, Acuity.max), fill = ANC.matched), offset = 0.005) + coord_fixed(ratio = 1) + theme_bw() + scale_fill_manual(values = (my.cols[(c(1,2,3,7,9))])) + ggtitle("Day 3") + annotate(geom="text", x=0.5, y=.8, label=paste0("τ =",round(cormatrix[rownames(cormatrix) == "ANC, D3",colnames(cormatrix) == "Acuity, Max"],3),"\npadj =",round(pmatrix[rownames(pmatrix) == "ANC, D3",colnames(pmatrix) == "Acuity, Max"],11)), color="black") + theme(legend.key.size=unit(4,'mm'), legend.text=element_text(size=5), legend.title=element_text(size=5), axis.title = element_text(size = 7))
p3 <- ggplot(data = metadata_temp[metadata_temp$Day == "D7",]) + geom_mosaic(aes(x = product(ANC.matched, Acuity.max), fill = ANC.matched), offset = 0.005) + coord_fixed(ratio = 1) + theme_bw() + scale_fill_manual(values = (my.cols[(c(1,2,3,7,9))])) + ggtitle("Day 7") + annotate(geom="text", x=0.5, y=.8, label=paste0("τ =",round(cormatrix[rownames(cormatrix) == "ANC, D7",colnames(cormatrix) == "Acuity, Max"],3),"\npadj =",round(pmatrix[rownames(pmatrix) == "ANC, D7",colnames(pmatrix) == "Acuity, Max"],11)), color="black") + theme(legend.key.size=unit(4,'mm'), legend.text=element_text(size=5), legend.title=element_text(size=5), axis.title = element_text(size = 7))

Figure 1C:

plot_grid(p1,p2,p3,ncol=3)

Now we begin to look at our transcriptomic data. We import the count and TPM data from RSEM, as well as the Ensembl ID to gene symbol conversions, and log-transform the TPM values:

Count <- read.xlsx(paste0(prefix,"Tables/TableS1.xlsx"), sheet = 7, rowNames = TRUE)
TPM <- read.xlsx(paste0(prefix,"Tables/TableS1.xlsx"), sheet = 8, rowNames = TRUE)

genepc <- read.delim(paste0(prefix,"Ensembl_to_Symbol.txt"))
logTPM <- log2(TPM + 1)
metadata_long <- merge(metadata_long, qc_data)
rownames(metadata_long) <- metadata_long$Public.Sample.ID

To get an overall sense of the data we start with a PCA of the protein-coding genes:

proteincoding <- genepc$Gene.stable.ID[genepc$Gene.type == "protein_coding"]
logPCG <- logTPM[which(rownames(logTPM) %in% proteincoding),]
logPCG <- logPCG[which(rowSums(logPCG) > 0),]
pcaResults <- prcomp(t(logPCG), center = T, scale. = T)
PCs <- data.frame(pcaResults$x[,1:2])
metadata_long$PC1 <- PCs$PC1
metadata_long$PC2 <- PCs$PC2

The main contributions to PC1 are Exonic rate and Median 3’ bias:

myPalette <- colorRampPalette((brewer.pal(9, "RdYlBu")))
sc <- scale_colour_gradientn(colours = rev(myPalette(100)), limits=c(0,70))
p1 <- ggplot(data = metadata_long, aes(x = -1*PC1, y = PC2, color = Exonic.Rate*100)) + geom_point(size = 0.9) + theme_bw() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), axis.text=element_text(size=12),axis.title=element_text(size=14), plot.title = element_text(lineheight=.8, face="bold", size = 16), axis.ticks.y=element_blank())  + ggtitle("PCA") + sc + coord_fixed(ratio = 2.2)
p1$labels$colour <- "Exonic Rate"

myPalette <- colorRampPalette((brewer.pal(9, "RdYlGn")))
sc <- scale_colour_gradientn(colours = rev(myPalette(100)), limits=c(min(metadata_long$Median.3..bias),1))
p2 <- ggplot(data = metadata_long, aes(x = -1*PC1, y = PC2, color = Median.3..bias)) + geom_point(size = 0.9) + theme_bw() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), axis.text=element_text(size=12),axis.title=element_text(size=14), plot.title = element_text(lineheight=.8, face="bold", size = 16), axis.ticks.y=element_blank())  + ggtitle("PCA") + sc + coord_fixed(ratio = 2.2)
p2$labels$colour <- "Median 3' Bias"

Figure S1B:

plot_grid(p1,p2)

We apply quality control filtration to the data and re-examine the PCA:

#Calculate mitochondrial percentages
mitorows <- grepl("^MT-",genepc$Gene.name)
mitoids <- genepc$Gene.stable.ID[mitorows]
mitocounts <- Count[rownames(Count) %in% mitoids,]
mitosums <- colSums(mitocounts)
totalcounts <- colSums(Count)
percent.mt <- mitosums/totalcounts*100
metadata_long$percent.mt <- percent.mt

#Perform quality control filtration based on RNASeQC parameters
metadata_filtered <- metadata_long[metadata_long$percent.mt < 20 & metadata_long$Genes.Detected > 10000 & metadata_long$Median.Exon.CV < 1 & metadata_long$Exon.CV.MAD < 0.75 & metadata_long$Exonic.Rate*100 > 25 & metadata_long$Median.3..bias < 0.9,]

logPCG_filtered <- logPCG[,colnames(logPCG) %in% rownames(metadata_filtered)]
logTPM_filtered <- logTPM[,colnames(logTPM) %in% rownames(metadata_filtered)]
TPM_filtered <- TPM[,colnames(TPM) %in% rownames(metadata_filtered)]
Count_filtered <- Count[,colnames(Count) %in% rownames(metadata_filtered)]

tf <- rowSums(TPM_filtered > 0.1) > ncol(TPM_filtered)*.2
TPM_filtered <- TPM_filtered[tf,]
Count_filtered <- Count_filtered[tf,]
logTPM_filtered <- logTPM_filtered[tf,]
tf <- rowSums(Count_filtered >= 6) > ncol(Count_filtered)*.2
TPM_filtered <- TPM_filtered[tf,]
Count_filtered <- Count_filtered[tf,]
logTPM_filtered <- logTPM_filtered[tf,]

#Repeat PCA
logPCG_filtered <- logPCG_filtered[which(rowSums(logPCG_filtered) > 0),]
pcaResults <- prcomp(t(logPCG_filtered), center = T, scale. = T)
PCs <- data.frame(pcaResults$x[,1:2])
metadata_filtered$PC1 <- PCs$PC1
metadata_filtered$PC2 <- PCs$PC2

myPalette <- colorRampPalette((brewer.pal(9, "RdYlBu")))
sc <- scale_colour_gradientn(colours = rev(myPalette(100)), limits=c(0,70))
p1 <- ggplot(data = metadata_filtered, aes(x = -1*PC1, y = -1*PC2, color = Exonic.Rate*100)) + geom_point(size = 0.9) + theme_bw() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), axis.text=element_text(size=12),axis.title=element_text(size=14), plot.title = element_text(lineheight=.8, face="bold", size = 16), axis.ticks.y=element_blank())  + ggtitle("PCA") + sc + coord_fixed(ratio = 1.3)
p1$labels$colour <- "Exonic Rate"

myPalette <- colorRampPalette((brewer.pal(9, "RdYlGn")))
sc <- scale_colour_gradientn(colours = rev(myPalette(100)), limits=c(min(metadata_filtered$Median.3..bias),1))
p2 <- ggplot(data = metadata_filtered, aes(x = -1*PC1, y = -1*PC2, color = Median.3..bias)) + geom_point(size = 0.9) + theme_bw() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), axis.text=element_text(size=12),axis.title=element_text(size=14), plot.title = element_text(lineheight=.8, face="bold", size = 16), axis.ticks.y=element_blank())  + ggtitle("PCA") + sc + coord_fixed(ratio = 1.3)
p2$labels$colour <- "Median 3' Bias"

Figure S1C:

plot_grid(p1,p2)

We check that quality control filtration has not introduced any biases in terms of patient demographics.

metadata_long$QC <- rownames(metadata_long) %in% rownames(metadata_filtered)
metadata_long$QC <- mapvalues(metadata_long$QC, from = c("FALSE","TRUE"), to = c("FAIL","PASS"))

histmaker <- function(parameter){
  metadata_temp <- metadata_long[,colnames(metadata_long) %in% c(parameter,"QC")]
  metadata_temp[,1] <- factor(metadata_temp[,1])
  histtable <- as.data.frame(matrix(0L, nrow = nlevels(metadata_temp[,1])*2, ncol = 3))
  colnames(histtable) <- c("Param","QC","Freq")
  histtable[,1] <- rep(levels(metadata_temp[,1]), each = 2)
  histtable[,2] <- rep(c("FAIL","PASS"),nlevels(metadata_temp[,1]))
  for (i in 1:nrow(histtable)){
     histtable$Freq[i] <- sum(metadata_temp[,1] == histtable$Param[i] & metadata_temp$QC == histtable$QC[i])
  }
  return(histtable)
}

parameter <- "COVID"
fishertest <- fisher.test(table(metadata_long[,colnames(metadata_long) == "QC"],metadata_long[,colnames(metadata_long) == parameter]))
histtable <- histmaker(parameter)
p1 <- ggplot(histtable, aes(x = Param, y = Freq, fill = QC)) + geom_bar(position = "dodge",stat = "identity") + theme_bw() + ylab("Count") + xlab(parameter) + coord_fixed(ratio = .0033) + scale_fill_manual(values = c("darkred","forestgreen")) + annotate("text", x=2, y=max(histtable$Freq), label= paste0("Fisher: ",round(fishertest$p.value, digits = 3)), colour = "black")

parameter <- "Acuity.max"
fishertest <- fisher.test(table(metadata_long[,colnames(metadata_long) == "QC"],metadata_long[,colnames(metadata_long) == parameter]))
histtable <- histmaker(parameter)
p2 <- ggplot(histtable, aes(x = Param, y = Freq, fill = QC)) + geom_bar(position = "dodge",stat = "identity") + theme_bw() + ylab("Count") + xlab(parameter) + coord_fixed(ratio = .018) + scale_fill_manual(values = c("darkred","forestgreen")) + annotate("text", x=2, y=max(histtable$Freq), label= paste0("Fisher: ",round(fishertest$p.value, digits = 3)), colour = "black")

parameter <- "severity.max"
fishertest <- fisher.test(table(metadata_long[,colnames(metadata_long) == "QC"],metadata_long[,colnames(metadata_long) == parameter]))
histtable <- histmaker(parameter)
p3 <- ggplot(histtable, aes(x = Param, y = Freq, fill = QC)) + geom_bar(position = "dodge",stat = "identity") + theme_bw() + ylab("Count") + xlab(parameter) + coord_fixed(ratio = .008) + scale_fill_manual(values = c("darkred","forestgreen")) + annotate("text", x=2, y=max(histtable$Freq), label= paste0("Fisher: ",round(fishertest$p.value, digits = 3)), colour = "black")

parameter <- "Age"
fishertest <- fisher.test(table(metadata_long[,colnames(metadata_long) == "QC"],metadata_long[,colnames(metadata_long) == parameter]))
histtable <- histmaker(parameter)
p4 <- ggplot(histtable, aes(x = Param, y = Freq, fill = QC)) + geom_bar(position = "dodge",stat = "identity") + theme_bw() + ylab("Count") + xlab(parameter) + coord_fixed(ratio = .024) + scale_fill_manual(values = c("darkred","forestgreen")) + annotate("text", x=2, y=max(histtable$Freq), label= paste0("Fisher: ",round(fishertest$p.value, digits = 3)), colour = "black")

#parameter <- "sex"
#fishertest <- fisher.test(table(metadata_long[,colnames(metadata_long) == "QC"],metadata_long[,colnames(metadata_long) == parameter]))
#histtable <- histmaker(parameter)
#p5 <- ggplot(histtable, aes(x = Param, y = Freq, fill = QC)) + geom_bar(position = "dodge",stat = "identity") + theme_bw() + ylab("Count") + xlab(parameter) + coord_fixed(ratio = .008) + scale_fill_manual(values = c("darkred","forestgreen")) + annotate("text", x=2, y=max(histtable$Freq), label= paste0("Fisher: ",round(fishertest$p.value, digits = 3)), colour = "black")
p5 <- ""

#parameter <- "race"
#fishertest <- fisher.test(table(metadata_long[,colnames(metadata_long) == "QC"],metadata_long[,colnames(metadata_long) == parameter]))
#histtable <- histmaker(parameter)
#p6 <- ggplot(histtable, aes(x = Param, y = Freq, fill = QC)) + geom_bar(position = "dodge",stat = "identity") + theme_bw() + ylab("Count") + xlab(parameter) + coord_fixed(ratio = .008) + scale_fill_manual(values = c("darkred","forestgreen")) + annotate("text", x=2, y=max(histtable$Freq), label= paste0("Fisher: ",round(fishertest$p.value, digits = 3)), colour = "black")
p6 <- ""

#parameter <- "ethnicity"
#fishertest <- fisher.test(table(metadata_long[,colnames(metadata_long) == "QC"],metadata_long[,colnames(metadata_long) == parameter]))
#histtable <- histmaker(parameter)
#p7 <- ggplot(histtable, aes(x = Param, y = Freq, fill = QC)) + geom_bar(position = "dodge",stat = "identity") + theme_bw() + ylab("Count") + xlab(parameter) + coord_fixed(ratio = .008) + scale_fill_manual(values = c("darkred","forestgreen")) + annotate("text", x=2, y=max(histtable$Freq), label= paste0("Fisher: ",round(fishertest$p.value, digits = 3)), colour = "black")
p7 <- ""

parameter <- "BMI"
fishertest <- fisher.test(table(metadata_long[,colnames(metadata_long) == "QC"],metadata_long[,colnames(metadata_long) == parameter]))
histtable <- histmaker(parameter)
p8 <- ggplot(histtable, aes(x = Param, y = Freq, fill = QC)) + geom_bar(position = "dodge",stat = "identity") + theme_bw() + ylab("Count") + xlab(parameter) + coord_fixed(ratio = .024) + scale_fill_manual(values = c("darkred","forestgreen")) + annotate("text", x=2, y=max(histtable$Freq), label= paste0("Fisher: ",round(fishertest$p.value, digits = 3)), colour = "black")

Figure S1E:

p1

p2

p3

p4

p8

We do observe that more samples are discarded from later time points:

parameter <- "Day"
fishertest <- fisher.test(table(metadata_long[,colnames(metadata_long) == "QC"],metadata_long[,colnames(metadata_long) == parameter]))
histtable <- histmaker(parameter)
p1 <- ggplot(histtable, aes(x = Param, y = Freq, fill = QC)) + geom_bar(position = "dodge",stat = "identity") + theme_bw() + ylab("Count") + xlab(parameter) + coord_fixed(ratio = .014) + scale_fill_manual(values = c("darkred","forestgreen")) + annotate("text", x=2, y=max(histtable$Freq), label= paste0("Fisher: p = ",round(fishertest$p.value, digits = 8)), colour = "black")

Figure S1F:

p1

This may be partially attributed to longer sample processing times.

Figure S1G:

wilcoxtest <- wilcox.test(Processing.Time ~ QC, metadata_long)
ggplot(metadata_long, aes(x = factor(QC), y = as.numeric(Processing.Time), fill = factor(QC))) + geom_boxplot(outlier.shape = NA) + geom_jitter() + theme_bw() + ylab("Processing Time (h)") + scale_fill_manual(values = c("darkred","forestgreen")) + coord_fixed(ratio = .08) + annotate("text", x=1, y=23, label= paste0("Wilcoxon: p = ",round(wilcoxtest$p.value, digits = 5)), colour = "black")

However, there are no significant differences in processing time between samples that pass or fail within the same day.

Figure S1H:

ggplot(metadata_long, aes(x = factor(Day), y = Processing.Time, fill = QC)) + geom_boxplot(outlier.shape = NA) + geom_point(position = position_jitterdodge(jitter.height = 0, jitter.width = 0.2), size = 0.5, alpha = .5) + theme_bw() + scale_fill_manual(values = c("darkred","forestgreen")) + xlab("Day") + ylab("Processing Time (h)") + coord_fixed(ratio = 0.21)

Next we want to get an estimate of the sample purity since we performed bulk RNA-seq of enriched neutrophils. In order to do so, we need a reference single-cell RNA-seq dataset for deconvolution using CIBERSORTx. We will use the Bonn Cohort 2 from Schulte-Schrepping et al. We start by importing the data and collapsing the cluster labels into major hematopoeitic cell lineages.

seuratwb <- readRDS(paste0(prefix,"seurat_COVID19_freshWB-PBMC_cohort2_rhapsody_jonas_FG_2020-08-18.rds"))

barcodecluster <- as.character(seuratwb$RNA_snn_res.0.8)
temp <- seuratwb$RNA_snn_res.0.8
names(barcodecluster) <- names(temp)
barcodecluster <- cbind(rownames(barcodecluster),barcodecluster)
barcodecluster <- as.data.frame(barcodecluster)
barcodecluster$barcode <- as.character(rownames(barcodecluster))
barcodecluster$barcode[barcodecluster$barcodecluster %in% c("0","1","4","8")] <- "Mature_Neutrophil"
barcodecluster$barcode[barcodecluster$barcodecluster %in% c("12","15")] <- "Immature_Neutrophil"
barcodecluster$barcode[barcodecluster$barcodecluster %in% c("3","10","14","18")] <- "Monocyte"
barcodecluster$barcode[barcodecluster$barcodecluster %in% c("2","5","6","13","16")] <- "T_NK"
barcodecluster$barcode[barcodecluster$barcodecluster %in% c("7","22")] <- "B"
barcodecluster$barcode[barcodecluster$barcodecluster %in% c("19")] <- "Plasmablast"
barcodecluster$barcode[barcodecluster$barcodecluster %in% c("9","11","17","20","21","23","24")] <- "Other"
seuratwb <- AddMetaData(seuratwb, metadata = subset(barcodecluster, select = c("barcode")), col.name = "Cell.Type")

Figure S2A:

DimPlot(seuratwb, reduction = "umap", group.by = "Cell.Type", label = FALSE) + scale_color_manual(values = c("grey60","tomato","forestgreen","skyblue","grey90","black","slateblue3")) + theme(axis.line.x = element_blank(), axis.line.y = element_blank(), axis.title.x = element_blank(), axis.ticks = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank())

For CIBERSORTx, we will generate a cell type signature matrix based on “pseudobulk” samples from a single cell type. We will add up the counts from all the cells of a given lineage for each patient and treat them as bulk samples.

pseudobulk <- matrix(0L, nrow = nrow(seuratwb@assays$RNA@counts), ncol = (length(levels(factor(seuratwb$donor))))*6)
rownames(pseudobulk) <- (seuratwb@assays$RNA@data@Dimnames[[1]])
idx <- (levels(factor(seuratwb$donor)))
pseudobulk <- as.data.frame(pseudobulk)
for (i in 0:(ncol(pseudobulk)/6 - 1)){
  k <- i*6+1
  colnames(pseudobulk)[k] <- paste(as.character(idx[i+1]),"_Mature_Neu",sep = "")
  colnames(pseudobulk)[k+1] <- paste(as.character(idx[i+1]),"_Immature_Neu",sep = "")
  colnames(pseudobulk)[k+2] <- paste(as.character(idx[i+1]),"_Mono",sep = "")
  colnames(pseudobulk)[k+3] <- paste(as.character(idx[i+1]),"_TNK",sep = "")
  colnames(pseudobulk)[k+4] <- paste(as.character(idx[i+1]),"_B",sep = "")
  colnames(pseudobulk)[k+5] <- paste(as.character(idx[i+1]),"_Plasmablast",sep = "")
}

#Patient BN-31 does not have all the cell types so they are excluded
for (i in 0:17){
  k <- i*6+1
  temp <- subset(seuratwb, subset = donor == as.character(idx[i+1]) & Cell.Type == "Mature_Neutrophil")
  pseudobulk[,k] <- rowSums(temp@assays$RNA@counts)
  temp <- subset(seuratwb, subset = donor == as.character(idx[i+1]) & Cell.Type == "Immature_Neutrophil")
  pseudobulk[,k+1] <- rowSums(temp@assays$RNA@counts)
  temp <- subset(seuratwb, subset = donor == as.character(idx[i+1]) & Cell.Type == "Monocyte")
  pseudobulk[,k+2] <- rowSums(temp@assays$RNA@counts)
  temp <- subset(seuratwb, subset = donor == as.character(idx[i+1]) & Cell.Type == "T_NK")
  pseudobulk[,k+3] <- rowSums(temp@assays$RNA@counts)
  temp <- subset(seuratwb, subset = donor == as.character(idx[i+1]) & Cell.Type == "B")
  pseudobulk[,k+4] <- rowSums(temp@assays$RNA@counts)
  temp <- subset(seuratwb, subset = donor == as.character(idx[i+1]) & Cell.Type == "Plasmablast")
  pseudobulk[,k+5] <- rowSums(temp@assays$RNA@counts)
}
for (i in 19:(ncol(pseudobulk)/6 - 1)){
  k <- i*6+1
  temp <- subset(seuratwb, subset = donor == as.character(idx[i+1]) & Cell.Type == "Mature_Neutrophil")
  pseudobulk[,k] <- rowSums(temp@assays$RNA@counts)
  temp <- subset(seuratwb, subset = donor == as.character(idx[i+1]) & Cell.Type == "Immature_Neutrophil")
  pseudobulk[,k+1] <- rowSums(temp@assays$RNA@counts)
  temp <- subset(seuratwb, subset = donor == as.character(idx[i+1]) & Cell.Type == "Monocyte")
  pseudobulk[,k+2] <- rowSums(temp@assays$RNA@counts)
  temp <- subset(seuratwb, subset = donor == as.character(idx[i+1]) & Cell.Type == "T_NK")
  pseudobulk[,k+3] <- rowSums(temp@assays$RNA@counts)
  temp <- subset(seuratwb, subset = donor == as.character(idx[i+1]) & Cell.Type == "B")
  pseudobulk[,k+4] <- rowSums(temp@assays$RNA@counts)
  temp <- subset(seuratwb, subset = donor == as.character(idx[i+1]) & Cell.Type == "Plasmablast")
  pseudobulk[,k+5] <- rowSums(temp@assays$RNA@counts)
}
pseudobulk <- pseudobulk[rowSums(pseudobulk) > 0,]
pseudobulk <- pseudobulk[-grep("\\.",rownames(pseudobulk)),]
pseudobulk <- pseudobulk[-grep("^MT",rownames(pseudobulk)),]

Then we convert the Ensembl gene IDs to symbols for the neutrophil data and filter the pseudobulk matrix to include only the overlapping genes between the two datasets:

genepctemp <- genepc[genepc$Gene.stable.ID %in% rownames(Count_filtered),]
Count_temp <- Count
for (i in 1:nrow(Count_temp)){
  if (length(which(genepctemp$Gene.stable.ID == rownames(Count_temp)[i])) == 1){
    if (!(genepctemp$Gene.name[which(genepctemp$Gene.stable.ID == rownames(Count_temp)[i])] %in% rownames(Count_temp))){
      rownames(Count_temp)[i] <- genepctemp$Gene.name[which(genepctemp$Gene.stable.ID == rownames(Count_temp)[i])]
    }
  }
}
pseudobulk_temp <- pseudobulk[rownames(pseudobulk) %in% rownames(Count_temp),]

The matrix that was used to generate the signature matrix is included in the Github.

Now we are able to run CIBERSORTx Generate Signature Matrix, setting limits of 50 to 100 genes per cell type and filtering for only hematopoietic genes. The CIBERSORTx cell type signature matrix:

Figure S2B:

matrix <- read.xlsx(paste0(prefix,"Tables/TableS1.xlsx"), sheet = 13, rowNames = TRUE)
par(mar=c(1,1,1,1))
heatmap3(matrix, Colv = NA, cexRow = 0.2)

Using this cell type signature matrix, we generate estimated cell type fractions in CIBERSORTx using default parameters. An overview of the cell type proportions:

rownames(genomic_signatures) <- genomic_signatures$Public.Sample.ID
cibersort_temp <- genomic_signatures[,colnames(genomic_signatures) %in% c("Mature_Neutrophil","Immature_Neutrophil","Monocyte","T_NK","B","Plasmablast","Neutrophil_total")]
cibersort_temp <- cibersort_temp[rownames(cibersort_temp) %in% rownames(metadata_filtered),]
cibersort_sorted <- cibersort_temp[order(cibersort_temp$Neutrophil_total),]
cibersort_sorted <- cibersort_sorted[,-7]
data_percentage <- t(cibersort_sorted*100)
coul <- c("forestgreen","tomato","skyblue","slateblue3","gray","black")

Figure S2C:

barplot(data_percentage, col=coul , border=NA, xlab=NA, axisnames = FALSE)

Search for differences in CIBERSORTx estimated cell fractions across COVID status for all samples:

metadata_filtered <- merge(metadata_filtered, genomic_signatures)
rownames(metadata_filtered) <- metadata_filtered$Public.Sample.ID
metadata_filtered$COVID <- mapvalues(metadata_filtered$COVID, from = c(0,1), to = c("Negative","Positive"))

my.cols <- brewer.pal(3, "Set2")

wilcoxtest <- wilcox.test(as.numeric(Neutrophil_total)*100 ~ factor(COVID), metadata_filtered)
p1 <- ggplot(metadata_filtered, aes(x = factor(COVID), y = as.numeric(Neutrophil_total)*100, fill = factor(COVID))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 3)), colour = "black") + ggtitle("Total")

wilcoxtest <- wilcox.test(as.numeric(Mature_Neutrophil)*100 ~ factor(COVID), metadata_filtered)
p2 <- ggplot(metadata_filtered, aes(x = factor(COVID), y = as.numeric(Mature_Neutrophil)*100, fill = factor(COVID))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 8)), colour = "black") + ggtitle("Mature")

wilcoxtest <- wilcox.test(as.numeric(Immature_Neutrophil)*100 ~ factor(COVID), metadata_filtered)
p3 <- ggplot(metadata_filtered, aes(x = factor(COVID), y = as.numeric(Immature_Neutrophil)*100, fill = factor(COVID))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 8)), colour = "black") + ggtitle("Immature")

wilcoxtest <- wilcox.test(as.numeric(Monocyte)*100 ~ factor(COVID), metadata_filtered)
p4 <- ggplot(metadata_filtered, aes(x = factor(COVID), y = as.numeric(Monocyte)*100, fill = factor(COVID))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 8)), colour = "black") + ggtitle("Monocyte")

wilcoxtest <- wilcox.test(as.numeric(T_NK)*100 ~ factor(COVID), metadata_filtered)
p5 <- ggplot(metadata_filtered, aes(x = factor(COVID), y = as.numeric(T_NK)*100, fill = factor(COVID))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 4)), colour = "black") + ggtitle("T/NK")

wilcoxtest <- wilcox.test(as.numeric(B)*100 ~ factor(COVID), metadata_filtered)
p6 <- ggplot(metadata_filtered, aes(x = factor(COVID), y = as.numeric(B)*100, fill = factor(COVID))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 3)), colour = "black") + ggtitle("B")

wilcoxtest <- wilcox.test(as.numeric(Plasmablast)*100 ~ factor(COVID), metadata_filtered)
p7 <- ggplot(metadata_filtered, aes(x = factor(COVID), y = as.numeric(Plasmablast)*100, fill = factor(COVID))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 10)), colour = "black") + ggtitle("Plasmablast")

Figure S2D:

plot_grid(p1,p2,p3,p4,p5,p6,p7,ncol = 7)

Search for differences in CIBERSORTx estimated cell fractions across disease severity for all samples regardless of COVID status:

metadata_filtered <- merge(metadata_filtered, genomic_signatures)
metadata_temp <- metadata_filtered[metadata_filtered$severity.max %in% c("severe","non-severe"),]

my.cols <- brewer.pal(3, "RdBu")

wilcoxtest <- wilcox.test(as.numeric(Neutrophil_total)*100 ~ factor(severity.max), metadata_temp)
p1 <- ggplot(metadata_temp, aes(x = factor(severity.max), y = as.numeric(Neutrophil_total)*100, fill = factor(severity.max))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[3],my.cols[1])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 18)), colour = "black") + ggtitle("Total")

wilcoxtest <- wilcox.test(as.numeric(Mature_Neutrophil)*100 ~ factor(severity.max), metadata_temp)
p2 <- ggplot(metadata_temp, aes(x = factor(severity.max), y = as.numeric(Mature_Neutrophil)*100, fill = factor(severity.max))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[3],my.cols[1])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 3)), colour = "black") + ggtitle("Mature")

wilcoxtest <- wilcox.test(as.numeric(Immature_Neutrophil)*100 ~ factor(severity.max), metadata_temp)
p3 <- ggplot(metadata_temp, aes(x = factor(severity.max), y = as.numeric(Immature_Neutrophil)*100, fill = factor(severity.max))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[3],my.cols[1])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 8)), colour = "black") + ggtitle("Immature")

wilcoxtest <- wilcox.test(as.numeric(Monocyte)*100 ~ factor(severity.max), metadata_temp)
p4 <- ggplot(metadata_temp, aes(x = factor(severity.max), y = as.numeric(Monocyte)*100, fill = factor(severity.max))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[3],my.cols[1])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 5)), colour = "black") + ggtitle("Monocyte")

wilcoxtest <- wilcox.test(as.numeric(T_NK)*100 ~ factor(severity.max), metadata_temp)
p5 <- ggplot(metadata_temp, aes(x = factor(severity.max), y = as.numeric(T_NK)*100, fill = factor(severity.max))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[3],my.cols[1])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 24)), colour = "black") + ggtitle("T/NK")

wilcoxtest <- wilcox.test(as.numeric(B)*100 ~ factor(severity.max), metadata_temp)
p6 <- ggplot(metadata_temp, aes(x = factor(severity.max), y = as.numeric(B)*100, fill = factor(severity.max))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[3],my.cols[1])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 3)), colour = "black") + ggtitle("B")

wilcoxtest <- wilcox.test(as.numeric(Plasmablast)*100 ~ factor(severity.max), metadata_temp)
p7 <- ggplot(metadata_temp, aes(x = factor(severity.max), y = as.numeric(Plasmablast)*100, fill = factor(severity.max))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[3],my.cols[1])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 4)), colour = "black") + ggtitle("Plasmablast")

Figure S2E:

plot_grid(p1,p2,p3,p4,p5,p6,p7,ncol = 7)

We repeat this analysis on Day 0 samples only as there were no COVID-19-negative samples collected after Day 0:

metadata_temp <- metadata_filtered[metadata_filtered$Day == "D0",]

my.cols <- brewer.pal(3, "Set2")

wilcoxtest <- wilcox.test(as.numeric(Neutrophil_total)*100 ~ factor(COVID), metadata_temp)
p1 <- ggplot(metadata_temp, aes(x = factor(COVID), y = as.numeric(Neutrophil_total)*100, fill = factor(COVID))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 3)), colour = "black") + ggtitle("Total")

wilcoxtest <- wilcox.test(as.numeric(Mature_Neutrophil)*100 ~ factor(COVID), metadata_temp)
p2 <- ggplot(metadata_temp, aes(x = factor(COVID), y = as.numeric(Mature_Neutrophil)*100, fill = factor(COVID))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 8)), colour = "black") + ggtitle("Mature")

wilcoxtest <- wilcox.test(as.numeric(Immature_Neutrophil)*100 ~ factor(COVID), metadata_temp)
p3 <- ggplot(metadata_temp, aes(x = factor(COVID), y = as.numeric(Immature_Neutrophil)*100, fill = factor(COVID))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 8)), colour = "black") + ggtitle("Immature")

wilcoxtest <- wilcox.test(as.numeric(Monocyte)*100 ~ factor(COVID), metadata_temp)
p4 <- ggplot(metadata_temp, aes(x = factor(COVID), y = as.numeric(Monocyte)*100, fill = factor(COVID))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 8)), colour = "black") + ggtitle("Monocyte")

wilcoxtest <- wilcox.test(as.numeric(T_NK)*100 ~ factor(COVID), metadata_temp)
p5 <- ggplot(metadata_temp, aes(x = factor(COVID), y = as.numeric(T_NK)*100, fill = factor(COVID))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 4)), colour = "black") + ggtitle("T/NK")

wilcoxtest <- wilcox.test(as.numeric(B)*100 ~ factor(COVID), metadata_temp)
p6 <- ggplot(metadata_temp, aes(x = factor(COVID), y = as.numeric(B)*100, fill = factor(COVID))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 3)), colour = "black") + ggtitle("B")

wilcoxtest <- wilcox.test(as.numeric(Plasmablast)*100 ~ factor(COVID), metadata_temp)
p7 <- ggplot(metadata_temp, aes(x = factor(COVID), y = as.numeric(Plasmablast)*100, fill = factor(COVID))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(wilcoxtest$p.value, digits = 15)), colour = "black") + ggtitle("Plasmablast")

Figure S3A:

plot_grid(p1,p2,p3,p4,p5,p6,p7,ncol = 7)

Next we want to understand determinants of sample purity. We can observe correlations between ANC and CIBERSORTx Total neutrophil fraction, as well as ALC and CIBERSORTx T/NK fraction.

metadata_temp <- metadata_filtered[metadata_filtered$COVID == "Positive",]
metadata_temp <- metadata_temp[complete.cases(metadata_temp$AMC.matched) & complete.cases(metadata_temp$ALC.matched),]

my.cols <- brewer.pal(3, "Set2")
p1 <- ggplot(metadata_temp, aes(x = factor(ANC.matched), y = (Neutrophil_total)*100)) + geom_boxplot(outlier.shape = NA, fill = "blue3") + geom_jitter(height = 0, alpha = 0.3) + theme_bw() + stat_compare_means() + ylab("CIBERSORTx Total Neutrophil (%)") + xlab("ANC Category") + scale_fill_manual(values = my.cols[c(1,2)]) + scale_y_continuous(limits = c(0,max(metadata_temp$Neutrophil_total)*100)) + coord_fixed(ratio = 1/14) + theme(legend.position = "none")

p2 <- ggplot(metadata_temp, aes(x = factor(ALC.matched), y = (T_NK)*100)) + geom_boxplot(outlier.shape = NA, fill = "mediumpurple4") + geom_jitter(height = 0, alpha = 0.3) + theme_bw() + stat_compare_means() + ylab("CIBERSORTx T/NK (%)") + xlab("ALC Category") + scale_fill_manual(values = my.cols[c(1,2)]) + scale_y_continuous(limits = c(0,max(metadata_temp$T_NK)*100)) + coord_fixed(ratio = 1/11) + theme(legend.position = "none")

p3 <- ggplot(metadata_temp, aes(x = factor(AMC.matched), y = (Monocyte)*100)) + geom_boxplot(outlier.shape = NA, fill = "skyblue") + geom_jitter(height = 0, alpha = 0.3) + theme_bw() + stat_compare_means() + ylab("CIBERSORTx Monocyte (%)") + xlab("AMC Category") + scale_fill_manual(values = my.cols[c(1,2)]) + scale_y_continuous(limits = c(0,max(metadata_temp$T_NK)*100)) + coord_fixed(ratio = 1/11) + theme(legend.position = "none")

Figure S3B:

plot_grid(p1,p2,p3,ncol=3)

We see that lower purity samples tend to have lower ANC and higher ALC. The converse is also true.

metadata_lo <- metadata_filtered[order(metadata_filtered$Neutrophil_total),]
metadata_lo <- metadata_lo[1:round(nrow(metadata_lo)*.25),]
comparisondf <- matrix(0L, nrow = 25, ncol = 3)
colnames(comparisondf) <- c("Var1","Var2","Freq")
comparisondf[,1] <- c(rep(5,5),rep(4,5),rep(3,5),rep(2,5),rep(1,5))
comparisondf[,2] <- c(rep(1:5,5))
for (i in 1:nrow(comparisondf)){
  comparisondf[i,3] <- sum(metadata_lo$ANC.matched == comparisondf[i,1] & metadata_lo$ALC.matched == comparisondf[i,2], na.rm = TRUE)
}
theme_nogrid <- function (base_size = 12, base_family = "") {
  theme_bw(base_size = base_size, base_family = base_family) %+replace% 
    theme(panel.grid = element_blank())   
}
my.cols <- brewer.pal(7, "Greens")
p1 <- ggplot(as.data.frame(comparisondf), aes(x = factor(Var1), y = Freq, fill = factor(Var2))) + geom_bar(position="stack", stat="identity") + scale_fill_manual(values = (my.cols[3:7])) + xlab("ANC Category") + ylab("Count") + theme_bw() + coord_fixed(ratio = .1)
p1$labels$fill <- "ALC Quintile"

metadata_hi <- metadata_filtered[order(metadata_filtered$Neutrophil_total),]
metadata_hi <- metadata_hi[round(nrow(metadata_hi)*.75):nrow(metadata_hi),]
comparisondf <- matrix(0L, nrow = 25, ncol = 3)
colnames(comparisondf) <- c("Var1","Var2","Freq")
comparisondf[,1] <- c(rep(5,5),rep(4,5),rep(3,5),rep(2,5),rep(1,5))
comparisondf[,2] <- c(rep(1:5,5))
for (i in 1:nrow(comparisondf)){
  comparisondf[i,3] <- sum(metadata_hi$ANC.matched == comparisondf[i,1] & metadata_hi$ALC.matched == comparisondf[i,2], na.rm = TRUE)
}
my.cols <- brewer.pal(7, "Greens")
p2 <- ggplot(as.data.frame(comparisondf), aes(x = factor(Var1), y = Freq, fill = factor(Var2))) + geom_bar(position="stack", stat="identity") + scale_fill_manual(values = (my.cols[3:7])) + xlab("ANC Category") + ylab("Count") + theme_bw() + coord_fixed(ratio = .1)
p2$labels$fill <- "ALC Quintile"

metadata_filtered$Purity <- -1*as.numeric(rownames(metadata_filtered) %in% rownames(metadata_lo)) + as.numeric(rownames(metadata_filtered) %in% rownames(metadata_hi))
metadata_filtered$Purity <- mapvalues(metadata_filtered$Purity, from = c(-1,0,1), to = c("lo","mid","hi"))
metadata_temp <- metadata_filtered[metadata_filtered$Purity %in% c("lo","hi"),]
fisher.test(table(metadata_temp$ANC.matched, metadata_temp$Purity))

## 
##  Fisher's Exact Test for Count Data
## 
## data:  table(metadata_temp$ANC.matched, metadata_temp$Purity)
## p-value < 2.2e-16
## alternative hypothesis: two.sided

fisher.test(table(metadata_temp$ALC.matched, metadata_temp$Purity))

## 
##  Fisher's Exact Test for Count Data
## 
## data:  table(metadata_temp$ALC.matched, metadata_temp$Purity)
## p-value = 0.01286
## alternative hypothesis: two.sided

Figure S3C:

p1

p2

Next, we view the CIBERSORTx estimated neutrophil fractions over time in COVID-19-positive patients:

my.cols <- brewer.pal(3, "Set2")
metadata_temp <- metadata_filtered[metadata_filtered$COVID == "Positive" & metadata_filtered$Day %in% c("D0","D3","D7"),]

kwtest <- kruskal.test(g = metadata_temp$Day, x = metadata_temp$Neutrophil_total)
p1 <- ggplot(metadata_temp, aes(x = factor(Day), y = as.numeric(Neutrophil_total)*100, fill = factor(Day))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2],my.cols[3])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(kwtest$p.value, digits = 15)), colour = "black") + ggtitle("Total")

kwtest <- kruskal.test(g = metadata_temp$Day, x = metadata_temp$Mature_Neutrophil)
p2 <- ggplot(metadata_temp, aes(x = factor(Day), y = as.numeric(Mature_Neutrophil)*100, fill = factor(Day))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2],my.cols[3])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(kwtest$p.value, digits = 37)), colour = "black") + ggtitle("Mature")

kwtest <- kruskal.test(g = metadata_temp$Day, x = metadata_temp$Immature_Neutrophil)
p3 <- ggplot(metadata_temp, aes(x = factor(Day), y = as.numeric(Immature_Neutrophil)*100, fill = factor(Day))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.08) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .04) + scale_fill_manual(values = c(my.cols[1],my.cols[2],my.cols[3])) + scale_y_continuous(limits = c(0,100), expand = c(0,0)) + annotate("text", x=1, y=95, label= paste0("p=",round(kwtest$p.value, digits = 34)), colour = "black") + ggtitle("Immature")

Figure 1D:

plot_grid(p1,p2,p3,ncol=3)

We also check if CIBERSORTx estimated neutrophil-related cell type fractions are associated with any clinical parameters.

my.cols <- brewer.pal(5,"Reds")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$Creatinine.matched),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p1 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(Creatinine.matched), y = Neutrophil_total, fill = factor(Creatinine.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p2 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(Creatinine.matched), y = Neutrophil_total, fill = factor(Creatinine.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p3 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(Creatinine.matched), y = Neutrophil_total, fill = factor(Creatinine.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)

my.cols <- brewer.pal(5,"Greens")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$CRP.matched),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p4 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(CRP.matched), y = Neutrophil_total, fill = factor(CRP.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p5 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(CRP.matched), y = Neutrophil_total, fill = factor(CRP.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p6 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(CRP.matched), y = Neutrophil_total, fill = factor(CRP.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)

my.cols <- brewer.pal(5,"Blues")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$Ddimer.matched),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p7 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(Ddimer.matched), y = Neutrophil_total, fill = factor(Ddimer.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p8 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(Ddimer.matched), y = Neutrophil_total, fill = factor(Ddimer.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p9 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(Ddimer.matched), y = Neutrophil_total, fill = factor(Ddimer.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)

my.cols <- brewer.pal(5,"Purples")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$LDH.matched),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p10 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(LDH.matched), y = Neutrophil_total, fill = factor(LDH.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p11 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(LDH.matched), y = Neutrophil_total, fill = factor(LDH.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p12 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(LDH.matched), y = Neutrophil_total, fill = factor(LDH.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)

# my.cols <- brewer.pal(10,"Paired")
# metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$intubated),]
# metadata_temp <- metadata_temp[metadata_temp$COVID == "1",]
# p13 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(intubated), y = Neutrophil_total, fill = factor(intubated))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = (my.cols[5:6]))
# p14 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(intubated), y = Neutrophil_total, fill = factor(intubated))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = (my.cols[5:6]))
# p15 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(intubated), y = Neutrophil_total, fill = factor(intubated))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = (my.cols[5:6]))

my.cols <- brewer.pal(10,"Paired")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$CXR.infiltrates),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p16 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(CXR.infiltrates), y = Neutrophil_total, fill = factor(CXR.infiltrates))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[7:8])
p17 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(CXR.infiltrates), y = Neutrophil_total, fill = factor(CXR.infiltrates))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[7:8])
p18 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(CXR.infiltrates), y = Neutrophil_total, fill = factor(CXR.infiltrates))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[7:8])

my.cols <- brewer.pal(12,"Paired")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$Trop_72h),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p19 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(Trop_72h), y = Neutrophil_total, fill = factor(Trop_72h))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[11:12])
p20 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(Trop_72h), y = Neutrophil_total, fill = factor(Trop_72h))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[11:12])
p21 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(Trop_72h), y = Neutrophil_total, fill = factor(Trop_72h))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[11:12])

Figure S4 (Total Neutrophil):

plot_grid(p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,ncol=3)

## Warning: Groups with fewer than two data points have been dropped.

## Warning: Groups with fewer than two data points have been dropped.

plot_grid(#p13,p14,p15,
          p16,p17,p18,p19,p20,p21,ncol=3)

my.cols <- brewer.pal(5,"Reds")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$Creatinine.matched),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p1 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(Creatinine.matched), y = Mature_Neutrophil, fill = factor(Creatinine.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p2 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(Creatinine.matched), y = Mature_Neutrophil, fill = factor(Creatinine.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p3 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(Creatinine.matched), y = Mature_Neutrophil, fill = factor(Creatinine.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)

my.cols <- brewer.pal(5,"Greens")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$CRP.matched),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p4 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(CRP.matched), y = Mature_Neutrophil, fill = factor(CRP.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p5 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(CRP.matched), y = Mature_Neutrophil, fill = factor(CRP.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p6 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(CRP.matched), y = Mature_Neutrophil, fill = factor(CRP.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)

my.cols <- brewer.pal(5,"Blues")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$Ddimer.matched),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p7 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(Ddimer.matched), y = Mature_Neutrophil, fill = factor(Ddimer.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p8 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(Ddimer.matched), y = Mature_Neutrophil, fill = factor(Ddimer.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p9 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(Ddimer.matched), y = Mature_Neutrophil, fill = factor(Ddimer.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)

my.cols <- brewer.pal(5,"Purples")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$LDH.matched),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p10 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(LDH.matched), y = Mature_Neutrophil, fill = factor(LDH.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p11 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(LDH.matched), y = Mature_Neutrophil, fill = factor(LDH.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p12 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(LDH.matched), y = Mature_Neutrophil, fill = factor(LDH.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)

# my.cols <- brewer.pal(10,"Paired")
# metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$intubated),]
# metadata_temp <- metadata_temp[metadata_temp$COVID == "1",]
# p13 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(intubated), y = Mature_Neutrophil, fill = factor(intubated))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = (my.cols[5:6]))
# p14 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(intubated), y = Mature_Neutrophil, fill = factor(intubated))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = (my.cols[5:6]))
# p15 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(intubated), y = Mature_Neutrophil, fill = factor(intubated))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = (my.cols[5:6]))

my.cols <- brewer.pal(10,"Paired")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$CXR.infiltrates),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p16 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(CXR.infiltrates), y = Mature_Neutrophil, fill = factor(CXR.infiltrates))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[7:8])
p17 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(CXR.infiltrates), y = Mature_Neutrophil, fill = factor(CXR.infiltrates))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[7:8])
p18 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(CXR.infiltrates), y = Mature_Neutrophil, fill = factor(CXR.infiltrates))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[7:8])

my.cols <- brewer.pal(12,"Paired")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$Trop_72h),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p19 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(Trop_72h), y = Mature_Neutrophil, fill = factor(Trop_72h))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[11:12])
p20 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(Trop_72h), y = Mature_Neutrophil, fill = factor(Trop_72h))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[11:12])
p21 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(Trop_72h), y = Mature_Neutrophil, fill = factor(Trop_72h))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[11:12])

Figure S4 (Mature Neutrophil):

plot_grid(p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,ncol=3)

## Warning: Groups with fewer than two data points have been dropped.

## Warning: Groups with fewer than two data points have been dropped.

plot_grid(#p13,p14,p15,
          p16,p17,p18,p19,p20,p21,ncol=3)

my.cols <- brewer.pal(5,"Reds")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$Creatinine.matched),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p1 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(Creatinine.matched), y = Immature_Neutrophil, fill = factor(Creatinine.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p2 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(Creatinine.matched), y = Immature_Neutrophil, fill = factor(Creatinine.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p3 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(Creatinine.matched), y = Immature_Neutrophil, fill = factor(Creatinine.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)

my.cols <- brewer.pal(5,"Greens")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$CRP.matched),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p4 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(CRP.matched), y = Immature_Neutrophil, fill = factor(CRP.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p5 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(CRP.matched), y = Immature_Neutrophil, fill = factor(CRP.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p6 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(CRP.matched), y = Immature_Neutrophil, fill = factor(CRP.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)

my.cols <- brewer.pal(5,"Blues")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$Ddimer.matched),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p7 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(Ddimer.matched), y = Immature_Neutrophil, fill = factor(Ddimer.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p8 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(Ddimer.matched), y = Immature_Neutrophil, fill = factor(Ddimer.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p9 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(Ddimer.matched), y = Immature_Neutrophil, fill = factor(Ddimer.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)

my.cols <- brewer.pal(5,"Purples")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$LDH.matched),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p10 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(LDH.matched), y = Immature_Neutrophil, fill = factor(LDH.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p11 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(LDH.matched), y = Immature_Neutrophil, fill = factor(LDH.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)
p12 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(LDH.matched), y = Immature_Neutrophil, fill = factor(LDH.matched))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 3) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols)

# my.cols <- brewer.pal(10,"Paired")
# metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$intubated),]
# metadata_temp <- metadata_temp[metadata_temp$COVID == "1",]
# p13 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(intubated), y = Immature_Neutrophil, fill = factor(intubated))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = (my.cols[5:6]))
# p14 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(intubated), y = Immature_Neutrophil, fill = factor(intubated))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = (my.cols[5:6]))
# p15 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(intubated), y = Immature_Neutrophil, fill = factor(intubated))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = (my.cols[5:6]))

my.cols <- brewer.pal(10,"Paired")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$CXR.infiltrates),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p16 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(CXR.infiltrates), y = Immature_Neutrophil, fill = factor(CXR.infiltrates))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[7:8])
p17 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(CXR.infiltrates), y = Immature_Neutrophil, fill = factor(CXR.infiltrates))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[7:8])
p18 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(CXR.infiltrates), y = Immature_Neutrophil, fill = factor(CXR.infiltrates))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[7:8])

my.cols <- brewer.pal(12,"Paired")
metadata_temp <- metadata_filtered[complete.cases(metadata_filtered$Trop_72h),]
metadata_temp <- metadata_temp[metadata_temp$COVID == "Positive",]
p19 <- ggplot(metadata_temp[metadata_temp$Day == "D0",], aes(x = factor(Trop_72h), y = Immature_Neutrophil, fill = factor(Trop_72h))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[11:12])
p20 <- ggplot(metadata_temp[metadata_temp$Day == "D3",], aes(x = factor(Trop_72h), y = Immature_Neutrophil, fill = factor(Trop_72h))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[11:12])
p21 <- ggplot(metadata_temp[metadata_temp$Day == "D7",], aes(x = factor(Trop_72h), y = Immature_Neutrophil, fill = factor(Trop_72h))) + geom_violin() + theme_bw() + scale_y_continuous(limits = c(0,1), expand = c(0,0)) + coord_fixed(ratio = 1.2) + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks = element_blank(), panel.grid = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + stat_compare_means(label.x = 1, label.y = .1) + scale_fill_manual(values = my.cols[11:12])

Figure S4 (Immature Neutrophil):

plot_grid(p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,ncol=3)

## Warning: Groups with fewer than two data points have been dropped.

## Warning: Groups with fewer than two data points have been dropped.

plot_grid(#p13,p14,p15,
          p16,p17,p18,p19,p20,p21,ncol=3)

We next generate UMAPs for sample visualization in Figure 1E. Figure S5 features projections of many clinical parameters and quality control parameters - as is the case with other parameters in the study, not all can be reproduced here.

set.seed(10101)
#tpm.umap <- umap(t(logTPM_filtered))
#metadata_filtered$umap1 <- tpm.umap$layout[,1]
#metadata_filtered$umap2 <- tpm.umap$layout[,2]

metadata_filtered$severity.covid <- matrix(0L, nrow = nrow(metadata_filtered), ncol = 1)
for (i in 1:nrow(metadata_filtered)){
  if (metadata_filtered$COVID[i] == "Positive"){
    metadata_filtered$severity.covid[i] <- metadata_filtered$severity.max[i]
  }
  if (metadata_filtered$COVID[i] == "Negative"){
    if (metadata_filtered$severity.max[i] == "H"){
      metadata_filtered$severity.covid[i] <- "H"
    }
    else {
      metadata_filtered$severity.covid[i] <- "non_COVID"
    }
  }
}

my.cols <- brewer.pal(3, "Set2")
p1 <- ggplot(metadata_filtered, aes(x = umap1, y = umap2)) + theme_bw() + geom_point(aes(colour = factor(COVID)),size = 3) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), panel.border = element_blank(), axis.title.x = element_blank(), axis.title.y = element_blank(), axis.ticks = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank()) + scale_color_manual(values = c(my.cols[1], my.cols[2]))
p1$labels$colour <- "COVID"

my.cols.2 <- brewer.pal(3, "RdBu")
p2 <- ggplot(metadata_filtered, aes(x = umap1, y = umap2)) + theme_bw() + geom_point(aes(colour = factor(severity.covid)),size = 3) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), panel.border = element_blank(), axis.title.x = element_blank(), axis.title.y = element_blank(), axis.ticks = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank()) + scale_color_manual(values = c("grey", my.cols[1], my.cols.2[3],  my.cols[2]))
p2$labels$colour <- "Severity"

a <- "Mature_Neutrophil"
myPalette <- colorRampPalette((brewer.pal(9, "RdYlBu")))
sc <- scale_colour_gradientn(colours = rev(myPalette(100)), limits=c(min(metadata_filtered[,colnames(metadata_filtered) == a])*100,max(metadata_filtered[,colnames(metadata_filtered) == a])*100)) 
p3 <- ggplot(metadata_filtered, aes(x = umap1, y = umap2, colour = (as.numeric(Mature_Neutrophil)*100))) + theme_bw() + geom_point(size = 3) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), panel.border = element_blank(), axis.title.x = element_blank(), axis.title.y = element_blank(), axis.ticks = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank()) + sc
p3$labels$colour <- "Mature Neu"

a <- "Immature_Neutrophil"
myPalette <- colorRampPalette((brewer.pal(9, "RdYlBu")))
sc <- scale_colour_gradientn(colours = rev(myPalette(100)), limits=c(min(metadata_filtered[,colnames(metadata_filtered) == a])*100,max(metadata_filtered[,colnames(metadata_filtered) == a])*100)) 
p4 <- ggplot(metadata_filtered, aes(x = umap1, y = umap2, colour = (as.numeric(Immature_Neutrophil)*100))) + theme_bw() + geom_point(size = 3) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), panel.border = element_blank(), axis.title.x = element_blank(), axis.title.y = element_blank(), axis.ticks = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank()) + sc
p4$labels$colour <- "Immature Neu"

Figure 1E:

plot_grid(p1,p2,p3,p4,ncol=2)

We perform differential expression to test for genes and pathways associated with COVID-19 infection. We use the helper script Neutrophil_DESeq2.R to import the function. As we will see, if we do not make a correction for plasmablast contamination, we will have many highly expressed immunoglobulin genes in the COVID+ side.

source(paste0(prefix,"Neutrophil_DESeq2.R"))

rownames(metadata_filtered) <- metadata_filtered$Public.Sample.ID
DESeq2_list <- Neutrophil_DESeq2(counts = Count_filtered, mdata = metadata_filtered, day = "D0")
dds <- DESeqDataSetFromMatrix(countData = DESeq2_list$Count_select, colData = DESeq2_list$coldata, design = ~ COVID)
dds <- DESeq(dds)

res <- as.data.frame(results(dds, name="COVID_Positive_vs_Negative"))
filenam <- "Day0_COVIDpos_vs_COVIDneg_uncorrected"
temp <- genepc[which(genepc$Gene.stable.ID %in% rownames(res)),]
res$symbol <- matrix(0L, nrow = nrow(res))
for (i in 1:nrow(res)){
  if (rownames(res)[i] %in% temp$Gene.stable.ID){
    res$symbol[i] <- temp$Gene.name[which(rownames(res)[i] == temp$Gene.stable.ID)]
  } else {
    res$symbol[i] <- rownames(res)[i]
  }
}
res$rank <- sign(res$log2FoldChange)*(-1)*log10(res$pvalue)
res <- res[complete.cases(res),]
res_sig <- res[res$padj < 0.05,]
#write.table(res,paste0(prefix,"DESeq2/",filenam,".txt"),sep = "\t")

resordered <- res[order(res$rank),]

log2fc <- as.numeric(resordered$log2FoldChange)
log10p <- as.numeric(-1*log10(resordered$pvalue))
padj <- resordered$padj
rank <- resordered$rank
symbol <- resordered$symbol
combo <- cbind(log2fc,log10p,padj,symbol,rank)
colnames(combo) <- c("log2fc","log10p","padj","symbol","rank")
rownames(combo) <- rownames(resordered)
combo <- as.data.frame(combo)
combo$rank <- as.numeric(combo$rank)
combo$log10p <- as.numeric(combo$log10p)
combo$log2fc <- as.numeric(combo$log2fc)
combo$significance <- as.numeric(combo$padj < 0.05)
combo$significance <- as.factor(combo$significance)

IGgenes <- c("IGKC","IGKV4-1","IGKV5-2","IGKV6-21","IGKV3D-20","IGKV3D-11","IGLV4-69","IGLV8-61","IGLV4-60","IGLV6-57","IGLV10-54","IGLV1-51","IGLV1-47","IGLV7-46","IGLV5-45","IGLV1-44","IGLV7-43","IGLV1-40","IGLV1-36","IGLV3-27","IGLV3-25","IGLV2-23","IGLV3-21","IGLV3-19","IGLV2-18","IGLV3-16","IGLV2-14","IGLV2-11","IGLV3-10","IGLV3-9","IGLV3-1","IGLC1","IGLC2","IGLC3","IGLC7","IGHA2","IGHG4","IGHG2","IGHA1","IGHG1","IGHG3","IGHD","IGHM","IGHV6-1","IGHV1-2","IGHV1-3","IGHV2-5","IGHV3-7","IGHV3-11","IGHV3-13","IGHV3-15","IGHV1-18","IGHV3-20","IGHV3-21","IGHV3-23","IGHV1-24","IGHV2-26","IGHV3-33","IGHV4-34","IGHV4-39","IGHV1-46","IGHV3-48","IGHV3-49","IGHV5-51","IGHV3-53","IGHV1-58","IGHV4-61","IGHV3-66","IGHV1-69","IGHV2-70D","IGHV3-73","IGLV9-49","IGHV3-64","IGKV3D-15","IGHV4-59","IGHV3-74","IGHV3-72","IGHV4-31","IGHV3-43","IGKV2D-30","IGKV1-6","IGKV3-20","IGKV1D-33","IGKV1-17","IGKV1-8","IGKV1-16","IGKV1D-16","IGKV2-24","IGKV3-11","IGKV1-9","IGKV1-33","IGKV1-39","IGKV2D-28","IGKV1D-17","IGKV2-30","IGKV2D-29","IGKV1-12","IGKV1-5","IGKV2-28","IGKV3-15","IGKV1-27","IGKV2D-40","IGKV1D-39","IGLVI-70","IGKV2-29","IGLL5","IGHV3-30","IGHV2-70","IGKV1D-13","IGHV4-4","IGLV2-8","IGHV1-69D","IGHV7-4-1","IGHV3-64D","IGHV5-10-1")
combo$color <- 0
combo$color[combo$symbol %in% IGgenes] <- 1

combo$labels <- 0
combo$labels[combo$symbol %in% c("MZB1","JCHAIN")] <- 1
combo$labels <- as.factor(combo$labels)

options(ggrepel.max.overlaps = Inf)
my.cols <- brewer.pal(3,"Set2")
plot1 <- ggplot(combo, aes(x = log2fc, y = log10p)) + geom_point(data = subset(combo, color == 1), colour = my.cols[2]) + geom_point(data = subset(combo, color == 0), colour = "grey") + geom_text_repel(data = subset(combo, labels == 1), aes(label = as.character(symbol)))  + theme_bw() + ylab("-Log10(p-value)") + xlab("Log2(FC)") + annotate("text", x=1.3, y=0, label= "COVID+", colour = my.cols[2]) + annotate("text", x=-1.2, y=0, label= "COVID-", colour = my.cols[1]) + coord_fixed(ratio = .13) + theme(panel.grid = element_blank())

Figure S6A:

plot1

Thus we will use CIBERSORTx estimated cell type fractions as well as an additional immunoglobulin score to regress non-neutrophil contamination out of the data.

source(paste0(prefix,"Pathway_scoring.R"))

IG.score <- Pathway_scoring(IGgenes)
metadata_filtered$IG.score <- (IG.score)

We check the correlations with this score with MZB1, a plasmablast marker, and MS4A1, a B cell marker.

goi <- "MZB1"
metadata_filtered$GOI <- t(logTPM_filtered[genepc$Gene.stable.ID[which(genepc$Gene.name == goi)],])
metadata_filtered$GOI <- as.numeric(metadata_filtered$GOI)
summry <- lm(GOI ~ IG.score, metadata_filtered)
p1 <- ggplot(metadata_filtered, aes(x = as.numeric(IG.score), y = GOI)) + geom_point() + theme_bw() + geom_smooth(method = "lm") + xlab("Immunoglobulin Score") + ylab("MZB1 Log2(TPM+1)") + annotate("text", x=-3, y=8.75, label= paste0("R2 =",round(summary(summry)$r.squared,3)), colour = "black") + annotate("text", x=-3, y=8, label= paste0("p =",round(summary(summry)$coefficients[,4],308)), colour = "black")

goi <- "MS4A1"
metadata_filtered$GOI <- t(logTPM_filtered[genepc$Gene.stable.ID[which(genepc$Gene.name == goi)],])
metadata_filtered$GOI <- as.numeric(metadata_filtered$GOI)
summry <- lm(GOI ~ IG.score, metadata_filtered)
p2 <- ggplot(metadata_filtered, aes(x = as.numeric(IG.score), y = GOI)) + geom_point() + theme_bw() + geom_smooth(method = "lm") + xlab("Immunoglobulin Score") + ylab("MS4A1 Log2(TPM+1)") + annotate("text", x=-3, y=5, label= paste0("R2 =",round(summary(summry)$r.squared,3)), colour = "black") + annotate("text", x=-3, y=4.5, label= paste0("p =",round(summary(summry)$coefficients[,4],18)), colour = "black")

Figure S6B:

p1

p2

View the immunoglobulin score on a UMAP:

a = "IG.score"
myPalette <- colorRampPalette((brewer.pal(9, "RdYlBu")))
sc <- scale_colour_gradientn(colours = rev(myPalette(100)), limits=c(min(metadata_filtered[,colnames(metadata_filtered) == a]),max(metadata_filtered[,colnames(metadata_filtered) == a]))) 
p1 <- ggplot(metadata_filtered, aes(x = umap1, y = umap2, colour = (as.numeric(IG.score)))) + theme_bw() + geom_point(size = 3) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), panel.border = element_blank(), axis.title.x = element_blank(), axis.title.y = element_blank(), axis.ticks = element_blank(), axis.text.x = element_blank(), axis.text.y = element_blank()) + sc
p1$labels$colour <- "IG Score"

Figure S6C:

p1

We check the IG score across COVID status, Day, and SeverityMax.

my.cols <- brewer.pal(3,"Set2")
p1 <- ggplot(metadata_filtered, aes(x = factor(COVID), y = IG.score, fill = factor(COVID))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.15) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .5) + scale_fill_manual(values = c(my.cols[1],my.cols[2])) + stat_compare_means()

Figure S6D:

p1

my.cols <- brewer.pal(5,"Set2")
p2 <- ggplot(metadata_filtered, aes(x = factor(Day), y = IG.score, fill = factor(Day))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.15) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .5) + scale_fill_manual(values = c(my.cols[1],my.cols[2],my.cols[3],my.cols[4],my.cols[5])) + stat_compare_means()

Figure S6E:

p2

my.cols <- brewer.pal(3,"RdBu")
p3 <- ggplot(metadata_filtered[metadata_filtered$severity.max %in% c("severe","non-severe"),], aes(x = factor(severity.max), y = IG.score, fill = factor(severity.max))) + geom_boxplot(outlier.shape = NA) + geom_jitter(alpha = 0.15) + theme_bw() + theme(legend.position = "none", axis.ticks.y = element_blank(), axis.text.y = element_blank(), plot.margin = unit(c(0, 0, 0, 0), "cm")) + xlab("") + ylab("") + coord_fixed(ratio = .5) + scale_fill_manual(values = c(my.cols[3],my.cols[1])) + stat_compare_means()

Figure S6F:

p3

Now we can perform the COVID+ vs. COVID- differential expression analysis.

DESeq2_list <- Neutrophil_DESeq2(counts = Count_filtered, mdata = metadata_filtered, day = "D0")
dds <- DESeqDataSetFromMatrix(countData = DESeq2_list$Count_select, colData = DESeq2_list$coldata, design = ~ Neutrophil_total + T_NK_factor + Monocyte_factor + IG_factor + Plasmablast_factor + COVID)
dds <- DESeq(dds)

res <- as.data.frame(results(dds, name="COVID_Positive_vs_Negative"))
filenam <- "Day0_COVIDpos_vs_COVIDneg_correct-NeuCont+TNK+Monocyte+Plasmablast+IG"
temp <- genepc[which(genepc$Gene.stable.ID %in% rownames(res)),]
res$symbol <- matrix(0L, nrow = nrow(res))
for (i in 1:nrow(res)){
  if (rownames(res)[i] %in% temp$Gene.stable.ID){
    res$symbol[i] <- temp$Gene.name[which(rownames(res)[i] == temp$Gene.stable.ID)] 
  } else {
    res$symbol[i] <- rownames(res)[i]
  }
}
res$rank <- sign(res$log2FoldChange)*(-1)*log10(res$pvalue)
res <- res[complete.cases(res),]
res_sig <- res[res$padj < 0.05,]
#write.table(res,paste0(prefix,"DESeq2/",filenam,".txt"),sep = "\t")

resordered <- res[order(res$rank),]

log2fc <- as.numeric(resordered$log2FoldChange)
log10p <- as.numeric(-1*log10(resordered$pvalue))
pvalue <- as.numeric(resordered$pvalue)
padj <- resordered$padj
rank <- resordered$rank
symbol <- resordered$symbol
combo <- cbind(log2fc,log10p,padj,pvalue,symbol,rank)
colnames(combo) <- c("log2fc","log10p","padj","pvalue","symbol","rank")
rownames(combo) <- rownames(resordered)
combo <- as.data.frame(combo)
combo$rank <- as.numeric(combo$rank)
combo$log10p <- as.numeric(combo$log10p)
combo$log2fc <- as.numeric(combo$log2fc)
combo$pvalue <- as.numeric(combo$pvalue)
combo$significance <- as.numeric(combo$padj < 0.05)
combo$significance <- as.factor(combo$significance)

combo$color <- 0
combo$color[combo$pvalue < 1e-4 & combo$log2fc > 0.5] <- 1
combo$color[combo$pvalue < 1e-4 & combo$log2fc < -0.5] <- -1

combo$labels <- 0
combo$labels[rownames(combo) %in% rownames(resordered)[1:20]] <- 1
combo$labels[rownames(combo) %in% rev(rownames(resordered))[1:20]] <- 1
combo$labels <- as.factor(combo$labels)

options(ggrepel.max.overlaps = Inf)
my.cols <- brewer.pal(3,"Set2")
plot1 <- ggplot(combo, aes(x = log2fc, y = log10p)) + geom_point(data = subset(combo, color == 0), colour = "grey") + geom_point(data = subset(combo, color == 1), colour = my.cols[2]) + geom_point(data = subset(combo, color == -1), colour = my.cols[1]) + geom_text_repel(data = subset(combo, labels == 1), aes(label = as.character(symbol)))  + theme_bw() + ylab("-Log10(p-value)") + xlab("Log2(FC)") + annotate("text", x=1.3, y=0, label= "COVID+", colour = my.cols[2]) + annotate("text", x=-1.2, y=0, label= "COVID-", colour = my.cols[1]) + coord_fixed(ratio = .21) + theme(panel.grid = element_blank())

Figure 1F:

plot1

Next we perform Gene Set Enrichment Analysis.

gmt.file <- gmtPathways(paste0(prefix,"all_gene_sets.gmt"))
res <- read.xlsx(paste0(prefix,"Tables/TableS1.xlsx"), sheet = 14)

ranking <- res[,"rank"]
names(ranking) <- res$symbol
set.seed(15001)
fgseaRes <- fgsea(pathways = gmt.file, 
                  stats = ranking,
                  minSize=25,
                  maxSize=1000,
                  eps = 0)
#write.table(fgseaRes[,1:7], file = paste0(prefix,"GSEA_",filenam,".txt"), sep = "\t")

The running enrichment score data points are generated using the following code.

getEnrichmentDataframe <- function (pathway, stats, gseaParam = 1, ticksSize = 0.2) {
    rnk <- rank(-stats)
    ord <- order(rnk)
    statsAdj <- stats[ord]
    statsAdj <- sign(statsAdj) * (abs(statsAdj)^gseaParam)
    statsAdj <- statsAdj/max(abs(statsAdj))
    pathway <- unname(as.vector(na.omit(match(pathway, names(statsAdj)))))
    pathway <- sort(pathway)
    gseaRes <- calcGseaStat(statsAdj, selectedStats = pathway, returnAllExtremes = TRUE)
    bottoms <- gseaRes$bottoms
    tops <- gseaRes$tops
    combo <- as.data.frame(cbind(tops,bottoms))
    combo$average <- matrix(0L, nrow = nrow(combo), ncol = 1)
    for (p in 1:nrow(combo)){
      combo$average[p] <- (combo$tops[p]+combo$bottoms[p])/2
    }
    combo <- cbind(combo,pathway)
    return(combo)
}

signalingpathways <- c("HALLMARK_INTERFERON_GAMMA_RESPONSE","HALLMARK_INTERFERON_ALPHA_RESPONSE","GO_CYTOKINE_MEDIATED_SIGNALING_PATHWAY","HALLMARK_TNFA_SIGNALING_VIA_NFKB","GO_INTERLEUKIN_1_BETA_PRODUCTION","GO_TOLL_LIKE_RECEPTOR_SIGNALING_PATHWAY","HALLMARK_IL6_JAK_STAT3_SIGNALING","GO_INTERLEUKIN_8_PRODUCTION","GO_TUMOR_NECROSIS_FACTOR_SUPERFAMILY_CYTOKINE_PRODUCTION","GO_INTERLEUKIN_6_PRODUCTION","GO_INSULIN_LIKE_GROWTH_FACTOR_RECEPTOR_SIGNALING_PATHWAY","GO_RESPONSE_TO_EPIDERMAL_GROWTH_FACTOR")

cellularpathways <- c("GO_DEFENSE_RESPONSE_TO_VIRUS","GO_INNATE_IMMUNE_RESPONSE","GO_REGULATION_OF_VIRAL_GENOME_REPLICATION","GO_OXIDATIVE_PHOSPHORYLATION","GO_HYDROGEN_PEROXIDE_CATABOLIC_PROCESS","GO_REGULATION_OF_CARBOHYDRATE_BIOSYNTHETIC_PROCESS","GO_RIBOSOME_ASSEMBLY","GO_COTRANSLATIONAL_PROTEIN_TARGETING_TO_MEMBRANE")

dataframe <- getEnrichmentDataframe(gmt.file[[signalingpathways[1]]], ranking)
signalingdf <- cbind(dataframe$pathway,dataframe$average,rep(signalingpathways[1],nrow(dataframe)))
colnames(signalingdf) <- c("rank","enrichment","pathway")
for (i in 2:length(signalingpathways)){
  dataframe <- getEnrichmentDataframe(gmt.file[[signalingpathways[i]]],
               ranking)
  temp <- cbind(dataframe$pathway,dataframe$average,rep(signalingpathways[i],nrow(dataframe)))
  colnames(temp) <- c("rank","enrichment","pathway")
  signalingdf <- rbind(signalingdf,temp)
}

dataframe <- getEnrichmentDataframe(gmt.file[[cellularpathways[1]]], ranking)
cellulardf <- cbind(dataframe$pathway,dataframe$average,rep(cellularpathways[1],nrow(dataframe)))
colnames(cellulardf) <- c("rank","enrichment","pathway")
for (i in 2:length(cellularpathways)){
  dataframe <- getEnrichmentDataframe(gmt.file[[cellularpathways[i]]],
               ranking)
  temp <- cbind(dataframe$pathway,dataframe$average,rep(cellularpathways[i],nrow(dataframe)))
  colnames(temp) <- c("rank","enrichment","pathway")
  cellulardf <- rbind(cellulardf,temp)
}

my.cols <- brewer.pal(12,"Paired")
p1 <- ggplot(as.data.frame(signalingdf), aes(x = as.numeric(rank), y = as.numeric(enrichment), colour = pathway)) + geom_point() + theme_bw() + scale_colour_manual(values = rev(my.cols)) + coord_fixed(ratio = 10000) + theme(panel.grid.minor = element_blank(), legend.text=element_text(size=6)) + xlab("Rank in Gene List") + ylab("Running Enrichment Score") + ggtitle("GSEA: Signaling Pathways")
my.cols <- brewer.pal(8,"Dark2")
p2 <- ggplot(as.data.frame(cellulardf), aes(x = as.numeric(rank), y = as.numeric(enrichment), colour = pathway)) + geom_point() + theme_bw() + scale_colour_manual(values = rev(my.cols)) + coord_fixed(ratio = 10000) + theme(panel.grid.minor = element_blank(), legend.text=element_text(size=6)) + xlab("Rank in Gene List") + ylab("Running Enrichment Score") + ggtitle("GSEA: Cellular Processes")

Figure 1G:

p1

Figure 1H:

p2

sessionInfo()

## R version 4.1.1 (2021-08-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] fgsea_1.18.0                umap_0.2.7.0               
##  [3] ggpubr_0.4.0                heatmap3_1.1.9             
##  [5] SeuratObject_4.0.2          Seurat_4.0.4               
##  [7] cowplot_1.1.1               ggmosaic_0.3.3             
##  [9] corrplot_0.90               openxlsx_4.2.4             
## [11] DESeq2_1.32.0               SummarizedExperiment_1.22.0
## [13] Biobase_2.52.0              MatrixGenerics_1.4.3       
## [15] matrixStats_0.60.1          GenomicRanges_1.44.0       
## [17] GenomeInfoDb_1.28.2         IRanges_2.26.0             
## [19] S4Vectors_0.30.0            BiocGenerics_0.38.0        
## [21] dplyr_1.0.7                 plyr_1.8.6                 
## [23] RColorBrewer_1.1-2          ggrepel_0.9.1              
## [25] ggplot2_3.3.5               knitr_1.33                 
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.2.2             reticulate_1.20        tidyselect_1.1.1      
##   [4] RSQLite_2.2.8          AnnotationDbi_1.54.1   htmlwidgets_1.5.3     
##   [7] grid_4.1.1             BiocParallel_1.26.2    Rtsne_0.15            
##  [10] munsell_0.5.0          codetools_0.2-18       ica_1.0-2             
##  [13] future_1.22.1          miniUI_0.1.1.1         withr_2.4.2           
##  [16] colorspace_2.0-2       highr_0.9              ROCR_1.0-11           
##  [19] ggsignif_0.6.2         tensor_1.5             listenv_0.8.0         
##  [22] labeling_0.4.2         GenomeInfoDbData_1.2.6 polyclip_1.10-0       
##  [25] bit64_4.0.5            farver_2.1.0           parallelly_1.27.0     
##  [28] vctrs_0.3.8            generics_0.1.0         xfun_0.25             
##  [31] fastcluster_1.2.3      R6_2.5.1               locfit_1.5-9.4        
##  [34] bitops_1.0-7           spatstat.utils_2.2-0   cachem_1.0.6          
##  [37] DelayedArray_0.18.0    assertthat_0.2.1       promises_1.2.0.1      
##  [40] scales_1.1.1           gtable_0.3.0           globals_0.14.0        
##  [43] goftest_1.2-2          rlang_0.4.11           genefilter_1.74.0     
##  [46] splines_4.1.1          rstatix_0.7.0          lazyeval_0.2.2        
##  [49] spatstat.geom_2.2-2    broom_0.7.9            yaml_2.2.1            
##  [52] reshape2_1.4.4         abind_1.4-5            backports_1.2.1       
##  [55] httpuv_1.6.2           tools_4.1.1            ellipsis_0.3.2        
##  [58] spatstat.core_2.3-0    ggridges_0.5.3         Rcpp_1.0.7            
##  [61] zlibbioc_1.38.0        purrr_0.3.4            RCurl_1.98-1.4        
##  [64] rpart_4.1-15           openssl_1.4.5          deldir_0.2-10         
##  [67] pbapply_1.4-3          zoo_1.8-9              haven_2.4.3           
##  [70] cluster_2.1.2          magrittr_2.0.1         data.table_1.14.0     
##  [73] RSpectra_0.16-0        scattermore_0.7        lmtest_0.9-38         
##  [76] RANN_2.6.1             fitdistrplus_1.1-5     hms_1.1.0             
##  [79] patchwork_1.1.1        mime_0.11              evaluate_0.14         
##  [82] xtable_1.8-4           XML_3.99-0.7           rio_0.5.27            
##  [85] readxl_1.3.1           gridExtra_2.3          compiler_4.1.1        
##  [88] tibble_3.1.4           KernSmooth_2.23-20     crayon_1.4.1          
##  [91] htmltools_0.5.2        mgcv_1.8-36            later_1.3.0           
##  [94] tidyr_1.1.3            geneplotter_1.70.0     DBI_1.1.1             
##  [97] MASS_7.3-54            Matrix_1.3-4           car_3.0-11            
## [100] igraph_1.2.6           forcats_0.5.1          pkgconfig_2.0.3       
## [103] foreign_0.8-81         plotly_4.9.4.1         spatstat.sparse_2.0-0 
## [106] annotate_1.70.0        XVector_0.32.0         stringr_1.4.0         
## [109] digest_0.6.27          sctransform_0.3.2      RcppAnnoy_0.0.19      
## [112] spatstat.data_2.1-0    Biostrings_2.60.2      rmarkdown_2.10        
## [115] cellranger_1.1.0       leiden_0.3.9           fastmatch_1.1-3       
## [118] uwot_0.1.10            curl_4.3.2             shiny_1.6.0           
## [121] lifecycle_1.0.0        nlme_3.1-152           jsonlite_1.7.2        
## [124] carData_3.0-4          viridisLite_0.4.0      askpass_1.1           
## [127] fansi_0.5.0            pillar_1.6.2           lattice_0.20-44       
## [130] KEGGREST_1.32.0        fastmap_1.1.0          httr_1.4.2            
## [133] survival_3.2-13        glue_1.4.2             zip_2.2.0             
## [136] png_0.1-7              bit_4.0.4              stringi_1.7.4         
## [139] blob_1.2.2             memoise_2.0.0          irlba_2.3.3           
## [142] future.apply_1.8.1